
# ============================================================================================================
# Replication Code for PRSM-Article "Strategic ambiguity of party positions in multiparty competition"
# by Thomas Braeuninger and Nathalie Giger, 02.02.2016
# ============================================================================================================

# note that is the replication code for the empirical analyses and the robustness test. The code for the 
# Monte Carlo simulation of the theoretical model (cp. Figure 1) is provided in the appendix of the paper


## clear memory
rm(list=ls())


## loading packages
require(foreign)
library(lme4)
library(texreg)
library(readstata13)
library(car)


## read data
d <- read.dta13("Braeuninger_Giger_2016.dta")


## coding
# party-system-year
country <- d$country
countrynew <- d$countrynew
countryname <- d$countryname
party <- d$party
year <- d$year
countryname[countrynew==211] <- "Belgium/Wallonia"
countryname[countrynew==212] <- "Belgium/Flanders"
countryname[countrynew==51] <- "UK"
countryshortname <- substr(countryname,1,3)
countryshortname[countrynew==211] <- "BeW"
countryshortname[countrynew==212] <- "BeF"
cy_name <- paste(countryshortname,year)
# cy = ID for party-system-year
cy <- as.integer(factor(countrynew*100000+year))

# platform mean: y1
# we use the wordscore estimate of the party position (using 2002 CMP data as reference scores)
y1 <- d$y4 
y1sq <- y1^2
# platform ambiguity: sd1
# we use the log of the standard deviation of the wordscore estimate
sd1 <- d$sd4
lnsd1 <- log(sd1)
# decisive party member: dpm
# we use the position of the mean party voter
dpm <- d$meanpartyvoter
# platform position and decisive party member: demeaned using mean voter of party-system-year
y11 <- y1 - d$meanvoter
dpmm <- dpm - d$meanvoter

# length of manifesto
N_all <- d$N_all
q <- log(1/N_all)

# vote share of parties (from CMP data)
pervote <- d$pervote

# niche party status ( a) party family, b) share left-right issues in manifesto)
nicheparty <- d$nicheparty
share <-d$niche_rile

# government participation
govpart <- d$govpart

# distance between mean voter and decisive party member
dist <- abs(d$meanvoter-d$meanpartyvoter)

# distance to closest competitor
min <-d$min

ddf = data.frame(cbind(
  sd1,lnsd1,y1,y1sq,countrynew,cy,q,pervote,nicheparty,govpart,share,dist,min,dpm))



#=====================================================
#Figure 2
#=====================================================

pdf(file="fig2.pdf", width = 12, height=9)
par(oma=c(0,0,0,0))
par(mar=c(6,6,6,5))
plot(y1,lnsd1, xlim=c(1,10), ylim=c(-5,-1),
     xlab="mean of party platform", 
     ylab="standard deviation of party platform (logged)", cex.lab=2,
     cex.axis=2, pch=21)
qq <- order(y1)
plx<-predict(loess(lnsd1 ~ y1, span=.7, degree=2), se=T)
lines(y1[qq],plx$fit[qq],lwd=2)
ci <- qt(0.975,plx$df)*plx$se
lines(y1[qq],plx$fit[qq] - ci[qq], lty=2,lwd=2)
lines(y1[qq],plx$fit[qq] + ci[qq], lty=2,lwd=2)
dev.off()


#=====================================================
#Table 1
#=====================================================

h1 <- lmer(y11 ~ dpmm + (1 + dpmm | cy), data=ddf)
summary(h1)
VarCorr(h1)
AIC(h1)
BIC(h1)
logLik(h1)
VarCorr(h1)
as.data.frame(VarCorr(h1))

texreg(h1, file="modelh1.tex", digits=3, stars=c(0.01, 0.05, 0.1), bookstabs=T,
       custom.coef.names = c("Intercept", "Decisive party member"))


#=====================================================
#Table 2
#=====================================================

# Model 1:
m1 <- lmer(lnsd1 ~ y1 + y1sq + q + (1 + y1 + y1sq | cy), 
           data=ddf, control = lmerControl(optimizer="Nelder_Mead"))
m1.robustness.check <- lmer(lnsd1 ~ y1 + y1sq + q + (1 + y1 + y1sq | cy),
           data=ddf, control = lmerControl(optimizer="bobyqa", optCtrl=list(maxfun=100000)))
summary(m1)
AIC(m1)
BIC(m1)
logLik(m1)
VarCorr(m1)
as.data.frame(VarCorr(m1))

# Model 2:
m2 <- lmer(lnsd1 ~ y1 + y1sq + q + pervote +  govpart + nicheparty + (1 + y1 + y1sq | cy),
           data=ddf, control = lmerControl(optimizer="Nelder_Mead"))
summary(m2)
AIC(m2)
BIC(m2)
logLik(m2)
VarCorr(m2)
as.data.frame(VarCorr(m2))

# Model 3:
m3 <- lmer(lnsd1 ~ y1 + y1sq + q + pervote +  govpart +  share + min + (1 + y1 + y1sq | cy),
           data=ddf, control = lmerControl(optimizer="Nelder_Mead"))
summary(m3)
AIC(m3)
BIC(m3)
logLik(m3)
VarCorr(m3)
as.data.frame(VarCorr(m3))

# table 2 output
texreg(list(m1,m2,m3), file="models123.tex", digits=3, stars=c(0.01, 0.05, 0.1), bookstabs=T,
       custom.coef.names = c("Intercept", "Platform mean", "Platform mean squared", "Wordcount",
                             "Party vote share","Government party", "Niche party", "Niche party (share RILE)",
                             "Distance to closest competitor"))


#=====================================================
# test U-shape relationship for model 1
#=====================================================

# interquartile range [r1,r2] = [4.728, 5.903]
# 2*r1 = 9.455, 2*r2 = 11.806
quantile(y1)
# We test for u-shape in the model 
#     lnsd1 = b0 + b1*y1 + b2*y1^2 + XB + epsilon
# at the limits of the IQR, r1 and r2
# the partial of lnsd1 w.r.t. y1 is b1 + 2*b2*y1
# Hypothesis then is H: b1 + 2*b2*r1 < 0 and b1 + 2*b2*r1 > 0
# The null Hypothesis is H0: b1 + 2*b2*r1 >= 0 or b1 + 2*b2*r1 <= 0
linearHypothesis(m1, "y1 + 9.455*y1sq = 0", verbose=TRUE)
linearHypothesis(m1, "y1 + 11.806*y1sq = 0", verbose=TRUE)
linearHypothesis(m1, c("y1 + 9.455*y1sq = 0", "y1 + 11.806*y1sq = 0"), verbose=TRUE)


#=====================================================
#Figure 3
#=====================================================

pdf(file="fig3.pdf", width=10, height=10)
par(oma=c(2,1,1,1)) 
par(mar=c(0,0,0,0))
par(mfrow=c(9,7))
cy_name_unique <- unique(cy_name)
x <- ddf$y1
a.hat.m1 <- coef(m1)$cy[,1]
b.hat.m1 <- coef(m1)$cy[,2]                
c.hat.m1 <- coef(m1)$cy[,3]                
d.hat.m1 <- coef(m1)$cy[,4]
for (j in 1:max(cy)) {
  plot (x[ddf$cy==j], ddf$lnsd1[ddf$cy==j],
        xlim=c(2.5,7.5), ylim=c(-5,-.7), xlab="", ylab="", cex.lab=1.2,
        cex.axis=1.1, pch=21, mgp=c(0,0,0), xaxt="n", yaxt="n", cex.main=1.1)
  text(3.9,-.9,cy_name_unique[j],cex=1.3, font=2)
  curve (a.hat.m1[j] + b.hat.m1[j]*x + c.hat.m1[j]*x*x + d.hat.m1[j]*mean(q[cy==j]), lwd=1.5, col="black", add=TRUE)
}
axis (1, c(2.5,7.5))
axis (4, c(-5,-1),las=1)
dev.off()


#=====================================================
# Robustness test: is variance not a function of the mean anyway?
#=====================================================

# run model with means and variance from a sample of 1000 fake texts.
# "samplemtexts.tab" provides point estimate and variance of wordscore 
# estimates for 1000 hypothetical texts
dat <- read.csv("sampledtexts.tab", sep=";")
sam_names <- dat[,1]
sam_y1 <- as.double(dat[,2])
sam_sd1<- log(sqrt(as.double(dat[,3])))

cnew <- countryname
cnew[countrynew==211] <- "Belgium_fr"
cnew[countrynew==212] <- "Belgium_fl"
cnew[countryname=="UK"] <- "United_Kingdom"
cnew[countryname=="Finland"] <- "Finland_fi"
cnew[countryname=="Luxembourg"] <- "Luxembourg_de"

# run a MC simulation N times
N <- 100
sim_sd1 <- rep(NA,321)
sim_y1 <- rep(NA,321)
sim_coef <- matrix(0,N,2)
sim_sign <- matrix(0,N,2)
sim_ushape <- matrix(0,N,2)
for (kk in 1:N) {
  for (j in 1:max(cy)) {
    n.parties <- length(cy[cy==j])
    case <- cnew[cy==j][1]
    sample <- sample(c(1:1000),n.parties)
    sim_sd1[cy==j] <- sam_sd1[sam_names==case][sample]
    sim_y1[cy==j] <- sam_y1[sam_names==case][sample]
  }
  sim_y1sq <- sim_y1^2
  mc.data <- as.data.frame(cbind(sim_sd1, sim_y1, sim_y1sq, cy))
  mc1 <- lmer(sim_sd1 ~ sim_y1 + sim_y1sq + (1 + sim_y1 + sim_y1sq | cy), data=mc.data,
              control = lmerControl(optimizer="bobyqa", optCtrl=list(maxfun=1000000)))
  tstat <- fixef(mc1) / sqrt(diag(vcov(mc1)))
  # save coefficients for mean and mean_sq as well as significance
  sim_coef[kk,1] <- fixef(mc1)[2]
  sim_coef[kk,2] <- fixef(mc1)[3]
  sim_sign[kk,1] <- dt(tstat[2],1)
  sim_sign[kk,2] <- dt(tstat[3],1)
  # save partial at IRQ (if u-shaped, these are neg & pos.)
  qq <- quantile(sim_y1,c(.05,.95))
  marg1 <- fixef(mc1)[2] + 2 * qq[1] * fixef(mc1)[3]
  marg2 <- fixef(mc1)[2] + 2 * qq[2] * fixef(mc1)[3]
  sim_ushape[kk,] <- cbind(marg1,marg2) 
}
# as NOT (marg1<0 & marg2>0), H0 can be rejected
plot(sim_ushape, xlim=c(-.1,.1), ylim=c(-.1,.1))
abline(h=0)
abline(v=0)



